Data from the test plate from IMR using “Universal Primers” from Parada et al, 2016. I followed the pipeline from Jesse McNichols and the Fuhrman lab protocols page.
Questions:
read=read.csv(file='/Users/oliviaahern/Documents/GitHub/TestPlate/raw_seqs.csv',header=T,row.names=1)
data=data.frame(read[,7:8])
#barplot(data$Raw)
How many reads in the community samples
sub=subset(read, Frac_Comm=="Comm")
par(mar=c(5,5,0,1))
barplot(sub$Raw,las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,log='x',
xlab="Log10 No. Raw Reads")
sub=subset(read, Frac_Comm=="Frac")
par(mar=c(5,5,0,1))
barplot(sub$Raw,las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,log='x',
xlab="Log No. Raw Reads")
plot(log10(sub$Raw), sub$Frac.RNA, xlab='Log10 No. of Raw Reads', ylab='cDNA concentration', pch =21, bg='black')
Percent abundance of prokaryotic, eukaryotic, and cyanobacterial reads after sorting raw sequences based on hits to Silva and PR2 databases
On average, 1.8% of the total sequences were Eukaryotes
On average 0.46% of the raw reads were eukaryotes, so about 163/47,000 reads.
reads=read.csv(file='/Users/oliviaahern/Documents/GitHub/TestPlate/pct_euk_bact.csv',header=T,row.names=1)
dim(reads)
## [1] 95 17
data1=data.frame(reads[,13:15])
sub=subset(reads, F_C=="Comm")
data1=data.frame(sub[,13:15])
par(mar=c(7,5,1,1),xpd=T)
barplot(as.matrix(t(data1)),las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,
xlab="% Abundance", col=c('navy','forestgreen','firebrick'),space=0)
box(which='plot')
legend(0,-10, legend=c("Prok", "Cyano", "Euks"), pch =22, pt.bg=c("firebrick", 'forestgreen','navy'),
bty='n', ncol=3)
On average, 3.4% of the fractionated sequences were eukaryotes, so about 806/24,000 reads.
Percentage of total raw reads that match to prokaryotic, eukarytoic, or cyanobacterial databases. The samples that don’t reach to 100 did not have matches to the databases, therefore we can assume they are probably chimera.s
sub=subset(reads, F_C=="Frac")
data1=data.frame(sub[,13:15])
barplot(as.matrix(t(data1)),las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,
xlab="% Abundance", col=c('navy','forestgreen','firebrick'),space=0)
box(which='plot')
legend(0,-10, legend=c("Prok", "Cyano", "Euks"), pch =22, pt.bg=c("firebrick", 'forestgreen','navy'),
bty='n', ncol=3)
Total of 110 unique ASVs.
library(phyloseq)
x<-read.csv(file='/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/feature-table.biom.csv',header=TRUE,row.names=1)
OTU = otu_table(x, taxa_are_rows=T)
taxa<-read.csv(file="/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/taxonomy.csv",header=TRUE,row.names=1)
t<-as.matrix(taxa)
tax2<-tax_table(t)
map<-import_qiime_sample_data("/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/map.txt")
#tree=read.tree("/Users/ahern/R/trophic_cascades/mc2/4Nov21/ASVs/tree.nwk")
phyo = phyloseq(OTU, tax2, map)
phyo
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 110 taxa and 84 samples ]
## sample_data() Sample Data: [ 84 samples by 10 sample variables ]
## tax_table() Taxonomy Table: [ 110 taxa by 9 taxonomic ranks ]
library(ggplot2)
#phyo_abund=subset_samples(phyo, Comm_Frac=="Comm")
phyo_abund=transform_sample_counts(phyo, function(x) x/sum(x))
pd <- psmelt(phyo_abund)
colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
#colors=readRDS('colors_class.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = ASV)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "none")
d_rgn
phy=subset_taxa(phyo, Class!="Fungi")
#phyo_abund=subset_samples(phyo, Comm_Frac=="Comm")
phyo_abund=transform_sample_counts(phy, function(x) x/sum(x))
pd <- psmelt(phyo_abund)
colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
#colors=readRDS('colors_class.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
## Warning: Removed 1680 rows containing missing values (position_stack).
#sub=subset_samples(phyo_abund, Comm_Frac=="Frac")
#write.csv(otu_table(sub),'18S_nofungi.csv')
#write.csv(tax_table(sub),'18S_nofungi_tax.csv')
sub=subset_samples(phyo_abund, Comm_Frac=="Frac")
control <- subset_samples(sub,SampleID == "Ace.TB5.1" | SampleID == "Ace.TB5.4" | SampleID == "Ace.TB5.5" | SampleID == "Ace.TB5.15" | SampleID == "C12.TB5.7" | SampleID == "C12.TB5.11" | SampleID == "Met.TB5.3" | SampleID == "Met.TB5.9" | SampleID == "Met.TB5.10" | SampleID == "Met.TB5.12")
pd <- psmelt(control)
colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
#colors=readRDS('colors_class.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, BD), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
library(dplyr);library(reshape2);library(ggplot2);library(cowplot);library(tidyverse)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ tibble 3.1.8 ✔ purrr 0.3.4
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
theme_set(theme_cowplot())
load("/users/oliviaahern/Downloads/Siders-tagseq-dataframe-26022020.RData",verbose=T)
## Loading objects:
## siders_tagseq_counts
head(siders_tagseq_counts)
## DATASET LOCATION SAMPLEID YEAR REP Feature.ID
## 87172 Axial ExtractControl CTRL 2019 2fa21026b4047e354f119a43c0341a65
## 87236 Axial ExtractControl CTRL 2019 6bbf315adb7368844fd1dc625be81e30
## 87273 Axial ExtractControl CTRL 2019 01d46e90107929f3903e4b5bff58a630
## 87616 Axial ExtractControl CTRL 2019 389950e760558e927385e40108737f9f
## 87957 Axial ExtractControl CTRL 2019 e613ce213c5added27988e9117e933fd
## 88999 Axial ExtractControl CTRL 2019 6d4984d3004dc93a767b711fce76c22a
## Taxon
## 87172 Eukaryota;Opisthokonta
## 87236 Eukaryota
## 87273 Eukaryota
## 87616 Eukaryota;Opisthokonta;Fungi;Basidiomycota;Agaricomycotina;Agaricomycetes;Sclerotium;Sclerotium_sp.;
## 87957 Eukaryota;Opisthokonta;Fungi;Basidiomycota;Pucciniomycotina;Cystobasidiomycetes;Rhodotorula
## 88999 Eukaryota
## Confidence level1 level2 level3 level4 level5
## 87172 0.7041164 Eukaryota Opisthokonta
## 87236 1.0000000 Eukaryota
## 87273 1.0000000 Eukaryota
## 87616 0.9742164 Eukaryota Opisthokonta Fungi Basidiomycota Agaricomycotina
## 87957 0.9971533 Eukaryota Opisthokonta Fungi Basidiomycota Pucciniomycotina
## 88999 1.0000000 Eukaryota
## level6 level7
## 87172
## 87236
## 87273
## 87616 Agaricomycetes Sclerotium;Sclerotium_sp.;
## 87957 Cystobasidiomycetes Rhodotorula
## 88999
## SAMPLE COUNT
## 87172 Axial_ExtractControl_CTRL_2019 1520
## 87236 Axial_ExtractControl_CTRL_2019 260
## 87273 Axial_ExtractControl_CTRL_2019 31
## 87616 Axial_ExtractControl_CTRL_2019 471
## 87957 Axial_ExtractControl_CTRL_2019 24
## 88999 Axial_ExtractControl_CTRL_2019 96
tax=tax_glom(phyo, taxrank="Phyla")
phyo_abund=transform_sample_counts(tax, function(x) x/sum(x))
control <- subset_samples(phyo,Timepoint=="0")
#write.csv(tax_table(control), 'controltax.csv')
#write.csv(otu_table(control), 'controlotu.csv')
pd <- psmelt(control)
#colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
colors <- c("#bfd3e6","#fa9fb5", "#74c476", "#fc8d59", "#807dba", "#238443","#e31a1c", "#ec7014", "#969696")
read=read.csv(file='controlotu.csv',header=T,row.names=1)
colors <- c("#bfd3e6","#fa9fb5", "#74c476", "#fc8d59", "#807dba", "#238443","#e31a1c", "#ec7014", "#969696")
dim(read)
## [1] 9 7
par(mar=c(15,5,1,10),xpd=T)
barplot(as.matrix(read), las=2, col = colors, space =0,
ylab="% Relative Abundance")
box(which='plot')
legend(8,.7, legend=row.names(read), pch=22,
pt.bg=colors)
Total of 1,632 ASVs.
library(phyloseq)
x<-read.csv(file='/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/16S/feature-table.biom.csv',header=TRUE,row.names=1)
OTU = otu_table(x, taxa_are_rows=T)
taxa<-read.csv(file="/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/16S/taxonomy_16s.csv",header=TRUE,row.names=1)
t<-as.matrix(taxa)
tax2<-tax_table(t)
map<-import_qiime_sample_data("/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/map.txt")
#tree=read.tree("/Users/ahern/R/trophic_cascades/mc2/4Nov21/ASVs/tree.nwk")
phyo = phyloseq(OTU, tax2, map)
phyo
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1632 taxa and 94 samples ]
## sample_data() Sample Data: [ 94 samples by 10 sample variables ]
## tax_table() Taxonomy Table: [ 1632 taxa by 8 taxonomic ranks ]
phyo1 = subset_taxa(phyo, !Order==" Chloroplast")
phyo2 = subset_taxa(phyo1, !Family==" Mitochondria")
phyo=phyo2
phyo2
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1608 taxa and 94 samples ]
## sample_data() Sample Data: [ 94 samples by 10 sample variables ]
## tax_table() Taxonomy Table: [ 1608 taxa by 8 taxonomic ranks ]
agg=tax_glom(phyo2, taxrank="Class")
t=transform_sample_counts(agg, function(x) x/sum(x))
sub=subset_samples(t, Comm_Frac=="Comm")
plot_bar(sub)
phyo_abund=subset_samples(t, Comm_Frac=="Comm")
pd <- psmelt(phyo_abund)
#colors=randomcoloR::randomColor(52)
#write_rds(colors, 'colors_class.rds')
colors=readRDS('/Users/oliviaahern/Documents/GitHub/TestPlate/colors_class.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
phyo_abund=subset_samples(t, Comm_Frac=="Frac")
pd <- psmelt(phyo_abund)
#colors=randomcoloR::randomColor(52)
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Frac), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning: Removed 104 rows containing missing values (position_stack).